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Energy dissipation and decoherence are at first glance harmful to acquiring long exciton lifetime desired for 
efficient photovoltaics. In the presence of both optically forbidden (namely, dark) and allowed (bright) excitons, 
however, they can be instrumental as suggested in photosynthesis. By simulating quantum dynamics of exciton 
relaxations, we show that the optimized decoherence that imposes a quantum-to-classical crossover with the 
dissipation realizes a dramatically longer lifetime. In an example of carbon nanotube, the exciton lifetime 
increases by nearly two orders of magnitude when the crossover triggers stable high population in the dark 
exciton. 


Sunlight is a clean, abundant, and sustainable energy 
source. Hence, the effective energy conversion from sunlight 
into electricity is a grand challenge in science and technol¬ 
ogy, which leads to an emerging interest on the post-silicon 
photovoltaic materials such as carbon nanotubes, quantum 
dots, transition metal dichalcogenides for their promising ap¬ 
plications CHI- The energy conversion consists of three 
processes of nonequilibrium quantum dynamics of excitons 
(bound electron-hole pairs): Exciton generation from a pho¬ 
ton, exciton energy transfer, and charge separation of the ex¬ 
citon into electrodes. In the first two processes, the efficient 
conversion requires both a high photon absorption rate and 
long exciton lifetime. However, the optimization is hampered 
by the reversibility between absorption and emission of a pho¬ 
ton in the quantum dynamics: While a high absorption rate 
of photon is desired in the exciton generation process, it also 
leads to a high charge recombination rate of exciton. 


Recently, photovoltaic models inspired by photosynthesis 
were studied from the viewpoint of stead y-st ate heat engines 
modified with discrete exciton states |ll6l[l7|]. A mechanism 
for the enhancement of the photocurrent was proposed by de¬ 
signing the ultrafast classical transition from bright to dark ex¬ 
citons within a hundred femtoseconds, where quantum coher¬ 
ence between photons and excitons and the resultant photolu¬ 
minescence causing the recombination were also ignored O. 
However, these classical descriptions of the photocell dynam¬ 
ics make the validity of the enhancement questionable because 
the inherent dichotomy between quantum coherence and de¬ 
coherence is a crucial issue in enhancing the quantum effi¬ 
ciency in the target time domain: On the one hand the de¬ 
sired ultrafast transition from the bright to the dark excitons 
necessarily requires quantum coherent dynamics, whereas on 
the other hand, the decoherence that imposes a quantum-to- 
classical crossover is assumed to immediately occur. 


In this work, we show that a desired remarkable enhance¬ 
ment of the exciton lifetime by simultaneously keeping high 
absorption rate is achieved by utilizing nonequilibrium energy 
dissipation and decoherence by phonons. Though the dissi¬ 
pation is in general an obstacle to high efficiency of energy 
conversion, there are cases where nonequilibrium dissipations 
and decoherence are actually instrumental in achieving long 
exciton lifetime because they can suppress the exciton recom¬ 
bination by making quantum dynamics irreversible via con¬ 
comitant quantum-to-classical crossover |0] . 


In fact, irreversible exciton dynamics is exploited in photo¬ 
synthesis that also includes the above mentioned three pro¬ 
cesses, namely, the exciton generation, the exciton energy 
transfer, and the charge separation iH H. Photosynthe¬ 
sis achieves a remarkably high quantum efficiency reaching 
nearly 100% ||3], which means that an absorbed photon is 
converted to an exciton with no recombination. In the ex¬ 
citon generation process, the absorbed energy is irreversibly 
transferred from the optically allowed bright exciton to an 
optically forbidden dark exciton Isi-fl^. which can act as a 
ratchet between the exciton generation process and the next 
energy transfer process where the quantum coherence plays a 
role again |[lll-[l5]. 


For the microscopic understanding of optimal energy con¬ 
versions, here we study the transient dynamics for the exciton 
generation process, taking account of the dichotomy in a uni¬ 
fied quantum manner. To this end, we construct and investi¬ 
gate an open quantum model that consists of bright and dark 
excitons coupled with phonons and a dissipative photon, with 
using the combined method of the generalized quantum mas¬ 
ter equation (GQME) ifTsIl and the quantum continuous mea¬ 
surement theory of photon counting (30. We found that 
realistic coexistence of coherence and decoherence is the key 
for high quantum efficiency of photovoltaics. The crossover 
from quantum to classical dynamics due to dissipation and de¬ 
coherence by phonons assists the long exciton lifetime: Rapid 
quantum energy transfer from the bright to dark excitons sup¬ 
presses the initial radiative loss, whereas the dark excitons be¬ 
come stable through the concomitant crossover. Our results 
indeed reveal why high efficiency of exciton generation from 
photons can be compatible with the low photoluminescence 
by violating the reversibility as observed experimentally. 

By taking a typical semiconducting single-walled carbon 
nanotube (SWCNT) with the chirality (6,5), as a model mate¬ 
rial, we demonstrate that the exciton lifetime becomes nearly 
two-orders of magnitude longer than the case without dark 
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FIG. 1: (Color online) Dissipative quantum model of excitons in 
semiconducting single-walled carbon nanotube (SWCNT) with the 
(6,5) chirality. Both the bright exciton generation and the charge re¬ 
combination result from the dipole coupling with photons with the 
same strength g. The photons are dissipated to environment with the 
rate of k. The bright exciton is coupled with two dark excitons via a 
phononcontinuum that causes energy dissipation and decoherence. 


excitons. Moreover, our results are consistent with the ex¬ 
perimental indications ll2ll - [^ . In the SWCNT, at least two 
dark exciton states exist with energy below the lowest bright 
one ll2ll - [^ : Even- parity d ark exciton and spin- 

triplet dark exciton Odllsdllssll . The weak couplings between 
the bright and the dark excitons can strongly affect the exciton 
dynamics after the photon absorption. Hence, we construct 
a Hamiltonian consisting of one bright and two dark exciton 
states coupled to a phononcontinuum and a single photon, as 
a minimal model of the SWCNT (see also Fig.[T]). The model 
Hamiltonian is given by 


iT = JTs + 


( 1 ) 


The first term is the Hamiltonian of the exciton-photon sys- 

with d\ and ^2 being the bosonic creation oper- 
ators for the photon, the bright exciton, and the two dark 
excitons (dark-1 and dark-2), respectively, where the dark- 
1 (dark-2) exciton has even-parity (triplet) symmetry. The 
bright exciton is coupled to the photon with the dipole cou¬ 
pling strength g. The second term represents the phononic 
bath, Hb = Z;r=dl,d 2 ^^rqplqprq whcrC pl^ TCads thc 
bosonic creation operator of a phonon with a momentum q. 
The energy of the photon, the excitons, and the phonons are 
represented by /icjph, ^r=di,d 2 , and respectively. The 

last term indicates the phonon-mediated coupling between the 
bright and dark excitons, ^int = Er=di,d 2 KrqiHKr + 
r){ptq + Prq) with the cffcctive coupling strength (. 

Simulating radiative lifetime of excitons requires non¬ 
unitary dissipations of the photon to environment or by photo 
detections. We simulate the non-unitary events by using quan¬ 
tum continuous measurement theo ry of p hoton counting with 


the photon dissipation rate, k >0 Ill9ll20|]. The photon count¬ 


ing gives the probability that m photons are absorbed to envi¬ 
ronments (or disappear at the instrument) until a time t > 0 
after the absorption of a photon (or the creation of the bright 
exciton) aXt = 0\ 


For feasible simulations of the excitons with the photon 
counting, we use the moment generating function (MGF) of 
the probability Ai{X; t) = X]m=-oo exp[imA]7^(m; t), with 
the conjugate variable A to the absorbed photon number m. 
The MGF dynamics can be explicitly determined by the fol¬ 
lowing equation of motion for a generalized density matrix 
with A: 


dtp{t; A) = —i[H^ p{t; X)]/h + Ke^^ap{t; A)a^ 

— K[a^ap{t; X)p{t; X)a^a], (2) 


where the MGF is given by the trace of p{t; A) as A4{X; t) = 
Tr[p(t; A)] for f > 0 with the initial boundary condition, p{t = 
0; A) = po. Here, po is the initial density matrix just after the 
photon absorption. Applying the Franck-Condon principle to 
the initial state at t = 0, the initial state is reasonably given 
by a separable state of the form po oc pg (8) exp[— 
where pg is the initial density matrix of the exciton-photon 
system with one bright exciton and T is the temperature of 
the phononic bath. 

From now on, we assume that the coupling strength Cdi,d2 
is weak and take the influence into account by the stan¬ 
dard second-order perturbation approximation in terms of 
^int ES]. In this approximation, the exciton dynamics is 
adequately described by the GQME for the generalized re¬ 
duced density matrix (RDM), p^(t; A) = TrB[p(t; A)] where 
TrB means the partial trace of the phononic bath, and the in¬ 
fluence of the phonons on the exciton system is represented 
by the spectral density Jr(cj) = — ^rq) for 

r = dl,d2. Here, we assume that Jr{oo) has the standard 
Ohmic form ifT^ . where 

7 ^ = 'sj^q ICrqp ^nd c^cut dcnotc the effective coupling 
strength between the exciton and phonons and the energy cut¬ 
off of the spectral density, respectively. 

Solving the GQME numerically, we calculate the general¬ 
ized RDM that gives the MGF as M{X] t) = Trs[p^(t; A)]. 
The probability is obtained from the inverse Fourier transform 
of the MGF. Note that the generalized RDM p^(t; A = 0) 
coincides with the actual density matrix of the reduced sys¬ 
tem where we ignore the counting history, p^{t) = p^(f; A = 
0) E IB a. Hereafter, we omit the script “S” repre¬ 
senting the subspace of the exciton-photon system such as 
p^(t) ^ p(t) for simplicity. Although the procedure is 
straightforward, the details of the formalism of the GQME 
are described in the Supplemental Material a. 

Actual time-resolved photoluminescence experiments have 
time-resolution limit At, which is introduced into our theory 
phenomenologically by the Gaussian averaging: V{m^t) = 
75=^ IZo '^) exp[-(T - i) V iV 2 AtZdT. From the 

time-averaged probability with At = 0.1 ps, the time- 
resolved photoluminescence, L{t), and the quantum yield at a 
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time t, Y(t), are given by 


m s 2: mdtV{m; t), 

m 

Y(<) ^ E mV{m] t)/mo. 


(3) 

(4) 


where mo = 1 is the initial number of the bright exciton at 
t = 0. Note that Y{t) represents the maximum quantum 
yield of photoluminescence experiments because the other 
environmental-photon-absorption is also included in Y (t). 

For a model calculation of (6,5) SWCNTs at room tem¬ 
perature, we choose the value of the parameters as follows: 
^<^ph = ^br = 1--27 eV, Sdi = 1.265 eV, ed 2 = 1-15 eV, 
T = 300 K, hujcut = 0.2 eV, Hg = 10.5 meV, H-fdi = 0.875 
meV, and h^d 2 = 0.25 meV. Effects of the environment are 
described by the single parameter, namely, the decay rate of 
photon, K,. Since k, strongly depends on the ambient solvents, 
matrices, and/or substrates, we vary k, over several orders of 
magnitude with keeping b: as a small parameter. The en¬ 
ergy levels of bright and dark excitons are estimated from the 
photoluminescence experiments 112 ll 13214351] . The cutoff fre¬ 
quency of the spectral density hujcut is determined from the 
density of states of phonons in SWCNTs ll3^ . The coupling 
strength 7 di,d 2 is determined from the numerical fitting that 
reproduces the experiment lil (see Fig. He)). The dipole 
coupling strength g is taken to be ^ 7di,d2- Note that one 

confirms that one order of magnitude difference in g does not 
affect our main results. 

First, we show the time-resolved photoluminescence L{t) 
for = 10 ns in Fig. Oa), which may be a typical result 
for SWCNTs in aqueous solutions. The simulated L{t) is ac¬ 
curately fitted by a tri-exponential function: 

Lm{t) = h exp{-t/Ti) + I2 exp{-t/T2) + hexpi-t/ra) 

(5) 

with Ti < r 2 < Ts. We note that similar multi -exp onential 
curves are also observed in the experiments ll2lil22ll4Qll4d]l . 
The decay constants obtained in our calculation differ by sev¬ 
eral orders of magnitude: The fast decay is characterized by 
ri = 65 ps while the intermediate and slow decay constants, 
r 2 and rs, are found at 890 ps and 1 gs, respectively. 

While the luminescence rapidly decreases within 1 ns as 
shown in Fig. Oa), the energy of the exciton-photon system 
E{t) = Tr[p{t)Hs] remains over 80% of sor even at 100 ns 
with the plateau from 5ns to 50 ns where the quantum yield 
holds lower than 10%. The significant difference in the de¬ 
cays of energy and luminescence qualitatively explains both 
the seemingly contradictory experiments of the time-resolved 
photoluminescence 122.1^ a nd the pump-probe transient 
absorption spectroscopy iili: Some time-resolved pho¬ 
toluminescence spectroscopies show the rapid luminescence 
decay with tens of ps and a low quantum yield lEniHlIi 
whereas a dark exciton survives over lOus in the pump-probe 
transient absorption spectroscopy 1 ^1^ ]. 

Exciton lifetime is a crucial factor characterizing the dy¬ 
namics. Here, we define the energy lifetime of excitons tlt 
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FIG. 2: (Color online) Transient exciton dynamics in a model of (6,5) 
SWCNT. (a) Time-resolved photoluminescence L{t) and quantum 
yield Y(t). The parameter of radiative dissipation to environment is 
fixed 2 Lt = 10 ns. The solid red line indicates the tri-exponential 
fitting of the numerical result, Lfit(/). The time-dependent energy 
E{t) is also plotted in the same figure, (b) Population dynam¬ 
ics of the exciton-photon system with nph(t) = Tr[p(t)Yd] and 
^r=br,di,d 2 (/) = Tr[p(t)N5r]. (c) Time evolution of von Neumann 
entropy S and the quantum mutual information /. The shadow region 
indicates the interval between the two peaks of the von Neumann en¬ 
tropy. 



as the energy decay time to 1 /e of the initial value for consis¬ 
tency with the single-exponential decay model: 

E{Ti^T)/E{0)=e-\ ( 6 ) 

From the definition, we obtain tlt = 880 ns, which is longer 
than forty times of that in the system with no dark exciton 
given by the inverse of the rate of the photon absorption to the 
environment = 2k,~^ = 20 ns as detailed later. Here, 

the factor “two” in 2k,~^ results from the halved residence 
time of the photon state by the Rabi oscillation between the 
photon and the bright exciton ll^ . 

The enhanced energy lifetime of the exciton-photon sys¬ 
tem is attributed to the fast irreversible relaxation pathway 
from the bright exciton to the dark excitons accompanied by 
quantum-to-classical crossover. We see the relaxation in the 
population dynamics in Fig. [2b). The initial dynamics has a 
quantum nature where the populations of the photon and the 
bright exciton are oscillating with the anti-phase with the fre¬ 
quency of g/2, which clearly indicates the Rabi oscillation. 
Then, the population is gradually transferred to the dark exci¬ 
tons after 10 ps with a reduction of the oscillation due to deco¬ 
herence effects by phonons. In particular, the dark-2 exciton 
is stabilized from a few ns to 100 ns with a high population 
via the quantum-to-classical crossover. 
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The irreversible crossover is confirmed by using the time- 
dependence of the von Neumann entropy and the quantum 
mutual information as shown in Fig.^c). The von Neumann 
entropy, S{t) = —Tr[p(t) Inp(f)], initially remains small, 
which means that the initial dynamics is dominated by the 
time-evolution of a quantum pure state. Then, the von Neu¬ 
mann entropy increases due to the decoherence by phonons, 
and shows two peaks approximately at 300 ps and 600 ns with 
the strong mixing of states. These mixings, however, have dif¬ 
ferent origins, i.e. the correlation between particles including 
quantum entanglement at ^ 300 ps and the classical stochas¬ 
tic mixing of the state at ^ 600 ns. 

To distinguish these two distinct origins of the mixing, we 
introduce quantum mutual information: The quantum mu¬ 
tual information provides a measure of correlation between 
subsystems of quantum states. Here, we calculate the quan¬ 
tum mutual information between the subsystem X consisting 
of the photon and the bright exciton and the subsystem Y 
of the dark excitons, I{t) = S{p{t)\\px{t) 0 pY{t))- Here 
S{p\\a) = Tr[p(ln(p) — ln((j))] is the quantum relative en¬ 
tropy that gives a measure of difference between two quan¬ 
tum states, and px{t) and pY{t) are the RDMs defined by 
the partial trace of the subsystem of Y and X, respectively: 
/3x(i) =TrY[/o(i)] andpY(i) =Trx[/o(i)]. 

While the correlation between the subsystems is enhanced 
approximately at the first peak of S{t), t :^300 ps, the quan¬ 
tum mutual information gives no peak at the second peak of 
S{t), which implies that the system becomes a nearly separa¬ 
ble state with essentially no quantum entanglement after the 
first peak as shown in Fig He). Hence, the crossover occurs 
near the first peak , which leads to the superselection of dark 
exciton states imposed by the decoherence [0]. Ultimately, 
the fast relaxation from the bright exciton to the dark excitons 
and the rapid quantum-to-classical crossover do not contradict 
each other in the time domain. 




T1 [ps] 


FIG. 3: (Color online) Environmental effects on photoluminescence, 
(a) Time-resolved photoluminescence with several choices of 
The solid red lines indicate the tri-exponential fitting functions on 
the numerical results, (b) A^-dependence of decay constants in tri¬ 
exponential fitting function, (c) Correlation between n and T 2 for 
several choices of , which ranges from 10 ps to 1 ps (black 
filled circle). The filled square (red) symbols indicate the exper¬ 
imentally measured decay constants for (6,5) SWCNT reproduced 
from Ref. IHj. 


Next, we examine the environmental effects on the time- 


resolved photoluminescence by monitoring k, as shown in 
Fig. Ha). The photoluminescence decays are well fit by a tri¬ 
exponential function irrespective of k, though the decay con¬ 
stants strongly depend on k,. The two decay constants ri and 
r 2 for the fast and intermediate decays, respectively, are satu¬ 
rated with an increase in as shown in Fig.Hb). Besides, 
the largest decay constant rs shows no saturation and propor¬ 
tional to 100/1:“^ with increasing These decay constants 
share the common origin: The positive correlation between 
ri and r 2 is clearly seen in Fig. He). The convex curve is 
consistent with the results of the experiment JH. 
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FIG. 4: (Color online) Environmental effects on energy relaxation. 
(a,b) Normalized energy Ejshr as a function of time and k, for the 
(6,5) SWCNT model (dissipative system), and for the hypothetical 
model with the same parameters as those in the dissipative system ex¬ 
cept for the energy levels, huph = £br = Sdi = Sd 2 (non-dissipative 
system). Solid lines indicate the contours of the normalized energy, 
(c) Energy lifetime tlt as a function of k in the dissipative and non- 
dissipative systems. The lifetime is normalized by 
that is the lifetime in the system with no dark exciton. (d) Normal¬ 
ized lifetime in the dissipative system as a function of the imaginary 
part of the dielectric function at ct; = Uph of the ambient dielectric 
medium, Im[e:(ct;ph)] = /^/cc;ph. 


Finally, we emphasize the importance of the energy dissipa¬ 
tion for the enhancement of the exciton lifetime. Figures H^) 
and (b) show the system energy E{t) as a function of time and 
K, for the systems where (a) the dark excitons have the lower 
energy compared to the bright one (dissipative system) or (b) 
all the excitons have the same energy level (non-dissipative 
system). The dissipative system obviously survives longer 
than the non-dissipative system for small k, where the irre¬ 
versible relaxation becomes dominant at the initial dynamics 
fort <100ps. 

The difference is quantified in the energy lifetime tlt 
as shown in Fig. He). Although the lifetime of the non- 
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dissipative system is enhanced compared to that in the system 
with no dark exciton (no-dark system) the 

enhancement factor is not large even for large The factor 

is about two. While, the lifetime in the dissipative system is 
about 45 times as long as for large k,~^. The large 

enhancement of the lifetime suddenly occurs within a small 
range of k. When all the energy loss to the environment can be 
represented by the dielectric loss of the ambient homogeneous 
linear dielectric medium, k is related with the background 
dielectric function e{uj) with uj being a frequency by using 
the Poyinting’s theorem in the classical electrodynamics S: 

= ^phlm[5(ccph)]. In the case, the sudden enhancement of 
lifetime is estimated to occur at Im[5(co’ph)] — 2 x 10“^ as 
shown in the inset of Fig. HJc) , which is comparable to that 
of water in room temperature, Im[5(co’ph)] — 1 x 10“^ ll^ . 
Therefore, the exciton lifetime of the SWCNT placed in the 
medium with smaller Im[5(u;ph)], e.g. air, can be strongly en¬ 
hanced. 

In summary, we have investigated the effect of dissipa¬ 
tion and decoherence on the long-term exciton dynamics in 
a model of the semiconducting (6,5) carbon nanotube. Coun¬ 
terintuitively, the dissipation and decoherence enhance the ex¬ 
citon lifetime in the relaxation pathway through the dark ex- 
citons, which paradoxically appears in experiments as a rapid 
decay of photoluminescence and a slow growth of quantum 
yield. In the search for efficient solar cells, such bad radia¬ 
tors are at first glance unpromising. However, we have shown 
this is not necessarily true. Our study confirms that semicon¬ 
ducting carbon nanotubes become potentially efficient pho¬ 
tovoltaic materials as suggested in the experiments ll3^ . and 
provides a further guideline and insight into high potentials of 
other unexplored materials. 
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SUPPLEMENTAL MATERIAL 

We describe the details of our method to calculate the den¬ 
sity matrix pit) and the probability V{m;t) referred to in the 
main article. Note that we use the notation of the doubled 
Hilbert space, i.e. the Liouville space, where a matrix a in the 
normal Hilbert space becomes a vector | a) with the Hilbert- 
Schmidt inner product, (pld") = Tr[p^a] 111 121] • A linear ma¬ 
trix operation becomes a linear map on the Liouville space, 
i.e. a superoperator. This description is useful for both the 
formulation and the numerical calculation in open quantum 
systems. 

Let us first consider the closed exciton system [Eq. (1) in 
the main article] with no interaction with the ambient environ¬ 
ment. In this case, the dynamics of the initial density matrix 
|po) is fully determined by the system Hamiltonian H. The 
time-evolution during the time interval [0, t] is given by the 
unitary superoperator, 

U{t) = exp[-i[£(^) - ^{H)]t/h], (SI) 

where £(•) (^(•)) is the map from any operator A in a Hilbert 
space to the superoperator £(A) (^(A)) in the Liouville space 
that makes the operator act to a matrix from left (right): 
£(A)|(j)) = I Ad)) and^(A)|d)) = |dA)). The density ma¬ 
trix after the time-evolution |p(t))) = f/(t)|po)) is consistent 
with the trivial results 

|/3(t))) = 2{exp[-iHt])^{exp[iHt])\po} 

= I exp[—exp['ii7t])). (S2) 

Note that in the first line in Eq. (IS^ . we use the commu¬ 
tative property between the superoperators mapped by £(•) 
and 2{A)^{B) = ^(5)£(A), and their homomorphic 
and antihomomorphic properties: A{AB) = 2{A)2{B) and 
£(A -\-B) = £(A) -h £(5), and ^{Ab) = ^{B)^{A) and 
^(A + 5) =^(A)+^(5). 

Erom now on, we shall consider the time-evolution in the 
open exciton system with photon dissipation to the environ¬ 
ment based on the continuous measurement theory of photon 
counting |3l0]. Assuming that the environment does not ab¬ 
sorb more than one photon during any infinitesimal time in¬ 
terval, we consider here only two fundamental processes: no¬ 
absorption and one-photon absorption. The one-photon ab¬ 
sorption process during any infinitesimal time interval 6t is 
represented by the superoperator J6t = hi6t2{a)^{aA with 
hi > 0 being the absorption rate in the presence of a free pho¬ 
ton, which converts an n-photon state to an (n — 1)-photon 
state iHQ]. Note that J6t is an adequate quantum dynamical 
map because it is linear and completely positive (CP) ffl. 
Since the normalized pre-absorption state is converted to 
the non-normalized post-absorption state by the immediate 
absorption, |p+)) = J(5t|p_)), the associated absorption prob¬ 
ability is given by Tr[p+] = ((/| J(5t|p_)) = hi6tTr[a^dp-]. 

Next, we consider the no-absorption process. When the 
superoperator Wq (t) represents the time-evolution under the 


condition that no photon is absorbed in the interval [0, f], the 
map lEo(-) should be continuous and homomorphic because 
the no-absorption process during [0, t] must be divided into a 
sequential process of several no-absorption processes: 

Wo{0)=I, (S3) 

Wo{ti^t2) = Wo{ti)Wo{t2), (S4) 

for alHi A 2 > 0, where I is the identity superoperator. These 
properties are the same as those of the unitary time-evolution 
IJ (t). The no-absorption process is, however, different from 
U(t) due to the presence of the environment. The canonical 
form is given by lii 

Wo{t) = exp[-i[£{H) - ^0)]t/h], (S5) 

% = H — thna) a/2. (S6) 

The anti-Hermitian part of H ensures a decrease in the asso¬ 
ciated no-absorption probability in the infinitesimal time in¬ 
terval 5t, (/|lEo((^t)|p_) = 1 — hi5tTx[d)dp-], which com¬ 
pensates for an increase in the photon-absorption probability 
= hi5tTr[d)dp-]. It is noteworthy that lEo(t) is 
the CP superoperator that satisfies Eq. ([S3]) and Eq. (IS4l) . and 
approaches U(t) in the limit of k, 0. 

The actual time-evolution superoperator W{t) is defined 
such that it represents a time-evolution of a density matrix 
during the time interval [0,t] with no information about the 
absorption. Accordingly, it is given by the sum of all possible 
combinations of the two fundamental processes: 

00 

W{t) = Wo{t)+J2w^it), (S7) 

m=l 



X Wo{t - tm)J " • l^o(^2 - (S8) 

where Wm{t) represents the m-photon absorption process 
during the time interval [0,t]- The actual time-evolution is 
also described by the recursive relation, 

W{t) = Wo{t)^ [ drW{t-r)JWo{r), (S9) 
Jo 

which clearly indicates that the actual time-evolution consists 
of the immediate one-photon absorption process and the no¬ 
absorption process. In addition, W{t) satisfies the following 
equation-of-motion (EOM): 

dtW{t) = [4 + >Cb + jCint]w{t). (SIO) 

where the Liouvillians are defined by jCs = —i[2{Hs) — 
^{Hs)]/h + /^£(a)^(at) - /^[£(ata) -f ^(ata)]/2, jCb = 
-i[£(^B) - ^{Hb)]/K and^Ant = -i[£(^int) - 

^{Hint)]/K respectively. Here ^s,B,int is defined in Eq. (1) 
of the main article. Accordingly, we can obtain the density 
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matrix after the actual time-evolution |p(t)} = W{t)\po'^ 
by solving the EOM. Note that the EOM has a Lindblad 
form isl] and |p(t))) is normalized as Tr[p{t)] = ((/|p(f))) = 
= Tr[po]. It is noted that if there is a photon 
decoupled from the excitons at f = 0 , the population decays 
exponentially with the lifetime of which is calculated 
analytically Jl]. 

Let us turn to the calculation of probability that m- 

photons are absorbed during the time interval [0, t]. Erom the 
above discussion, we have already obtained the superoperator 
that represents the m-photon absorption process during [ 0 , t] 

Vim^t) = {i\Wmit)\po}- (Sll) 


Because the time-evolution equation (IS13I) has a form very 
similar to that for the standard density matrix, we make use of 
an established method developed in the quantum master equa¬ 
tion (QME). Then, we derive a closed EOM for the general¬ 
ized reduced density matrix in the exciton-photon subspace 
\p^{t; A)))g defined by 

|;5S(f;A)))g = B((/B|/5(f;A))), (S14) 

where b((^b | means the partial trace of the phonon bath sub- 

space b((Jb Ip) = |TrB[/|/5]))g = |TrB[p]))g. In accordance 
with the standard second-order time-convolutionless projec¬ 
tion operator method we start with the density matrix in 
the interaction picture. 


Although it is difficult to calculate the probability analytically 
due to the phonon interaction, we are able to numerically cal¬ 
culate it with an approximation based on the moment generat¬ 
ing function (MGE) of the probability, 

oo 

Af (A; t)=Y, t) = A))), (S12) 

m=0 

where we define the generalized density matrix by | p{t; A))) = 
W{t-,^)\po} with = ^(f)lj^jexp(a) being the 

time-evolution superoperator. The EOM for |p(t; A))) is there¬ 
fore given by the following equation, 

= [>Cs(A) + >Cb + Ant] |p(t; A))), (S13) 

with jCs{X) = >^s| j^jexp(iA)- boundary condition is 
given by the actual density matrix att = 0: |p(0;A))) = 
|po}- Here, the initial density matrix |A)) is assumed to 
be the tensor product state of a bright exciton state and the 
phonon bath in the thermal equilibrium state, which rea¬ 
sonably describes the state just after the exciton generation 
in accordance with the Eranck-Condon principle: |po} = 

\Po}s ® |/5?q))B ^0 = ^L|0)(0|V and = 

ex.p[—H b / ksT] /TrB[ex.p[—H b / ksT]]. Here, TrB is the par¬ 
tial trace with respect to the subspace of the phonon bath, ks 
is the Boltzmann constant, and T is the temperature of the 
bath. 

Equation (IS13b is equivalent to Eq. (2) in the main arti¬ 
cle, which is written in the ordinary Hilbert space. Because 
the actual time-evolution of the density matrix is given by 
|p(t))) = |p(t;0))) [compare Eq. (ISlOb and Eq. (IS13b 1. we 
hereafter focus on the calculation of the EOM (IS13b . The 
probability V{m;t) is given by the inverse Eourier transfor¬ 
mation of the MGE. 

In the main article, we assume that the coupling between 
the phonons and the excitons is weak. Hence, we include the 
phonon effects on the exciton-photon system within the sec¬ 
ond order perturbation in terms of the interaction Liouvillian 
Cint . The phonon bath is also assumed to remain in the ther¬ 
mal equilibrium at all time. 


|pi(f; A))) = exp[-(£s(A) -h ^b)^] |p(f; A))) (S15) 

Because |pi(f; A)) satisfies the following EOM, 

dt\pi{t;X)} = 4nt(^;^)|pi(i;^)»> (S16) 

with £{„,(*; = exp[-(£s(A) -f £B)f]Ant(f) exp[(£s(A) -f 

Zb)^], the formal solution is given by 

|pi(i; A)» = f exp[ f A)dr]|po)). (S17) 

Jo 

Introducing the explicit coupling constant a to the Liouvillian, 
aC\^^{r; A), we expand the above formal solution within the 
second order |[§] . Then, we obtain the following EOM for the 
generalized reduced density matrix in the interaction picture 
|/5f(i;A)))g = B((4|pi(f; A))), 

dt\pf{t;X)}^ = exp[-£s(A)f]t(f; A) exp[£s(A)f]|pf(f; A)»g, 

(SI 8) 


where we define 


nt POO 

t(f;A) = -/ dr du) 

- - AO 


r=dl,d 2 


huj 


X [coth(^^j^)cos(wT)S^ (t;A) - i sin(wr)B+(T; A)]. 

(S19) 


Here, = £(S]; 6 br + ± + ^Ir^r) and 

B^{t; A) = exp[Cs{X)t]B^ exp[—£s(A)f]. The influence of 
phonons on the excitons is represented by the spectral density 
Jr{uj) = — ftrq) for r = dl, d 2 , which is here 

assumed to have an Ohmic form; (lo) = 2j^L0P(L0)P(L0cut- 
‘^)/‘^cut where we introduce a cutoff frequency Wcut and the 
effective coupling strength 7 ^. 

Converting the equation in the interaction picture to that in 
the Schrodinger picture, we obtain the generalized QME for 
numerical calculation: 

dtlP(t;X)}g = [4(A) + t(f;A)]|p®(f;A)))g. (S20) 
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The initial condition is given by |p^(0; A)))g = |po)s- 
using (t; A), the MGF and the actual reduced density matrix 
of the exciton-photon subsystem (t) are calculated as 

p^(t)=p^(t;0), (S21) 

M(A;t) = (/>(t;A)» = s((/s|p"(t;A)»g. (S22) 

We solve the generalized QME (IS20b numerically. In the 
Liouville space where a density matrix is vectorized and su¬ 
peroperators can be represented by matrices, the QME (IS20b 
becomes just a system of linear differential equations. Hence, 
we use ordinary methods: Fourth-order Runge-Kutta (RK4) 
method and eigen-decomposition (ED). At the initial time, 
T(t; A) strongly depends on the time t. Therefore, we solve 
the equation by the RK4 method until T(t; A) becomes time- 
independent. We take the unit time step typically as 0.005 ps. 
After sufficiently long time, T(t; A) hardly changes in time 
and becomes a nearly time-independent superoperator, which 
is known as the Markovian limit. In the present calculations, 
the convergent time tc is typically tens of ps and is shorter 
than 100 ps. After the convergence, we carry out the ED of 
the generator Cc{X) = jCs{X) + T(tc; A): 

C^{X) = P{X)D{X)p-\X) (S23) 

where P(A) is the matrix composed of the eigenvectors, and 


D{X) is the diagonal matrix constructed from the correspond¬ 
ing eigenvectors. The generalized reduced density matrix 
p^{t;X) at any time after tc is directly given by 

X)}=P{X)exp[i){X){t-tc)]p-^{X)\p{tc;X)}. 

(S24) 

The method combining RK4 and ED works efficiently for the 
long-term dynamics calculation over several orders of time. 
We have confirmed that the results calculated from the com¬ 
bined method are consistent with those obtained only from the 
RK4 method in the interval [0,1 ns]. 
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